Test- Bed Simulations of Collisionless, Self-Gravitating Systems 

Using the Schrodinger Method 



ON 
ON 



in 



George Davies and Lawrence M. Widrow 
Department of Physics, Queen's University, Kingston, K7L 3N6, Canada 

February 5, 2008 



Abstract 



The Schrodinger Method is a novel approach for modeling numerically self-gravitating, col- 
lisionless systems that may have certain advantages over N-body and phase space methods. In 
particular, smoothing is part of the dynamics and not just the force calculation. This paper 
describes test-bed simulations which illustrate the viability of the Schrodinger Method. We de- 
velop the techniques necessary to handle "hot" systems as well as spherically symmetric systems. 
We demonstrate that the method can maintain a stable, equilibrium star cluster by constructing 
and then evolving a Plummer sphere. We also consider nonequilibrium initial conditions and 
follow the system as it attempts to reach virial equilibrium through phase mixing and violent 
relaxation. Finally we make a few remarks concerning the dynamics of axions and other bosonic 
\Q , dark matter candidates. The Schrodinger Method, in principle, provides an exact treatment of 

ON ■ these fields. However such "scalar field" simulations are feasible and warrented only if the de- 

r-j . Broglie wavelength of the particle is comparable to the size of the system of interest, a situation 

O that is almost certainly not the case for axions in the Galaxy. Therefore the Schrodinger Method 

treats axions in the same way as other collisionless particles. We challenge recent claims in the 
■ literature that axions in the Galaxy form soliton stars. 

CO 

1 Introduction 

Recently Widrow and Kaiser (1993) have developed a new numerical method for modeling self- 
gravitating, collisionless systems. Known as the Schrodinger Method (SM) this approach allows 
one to follow the evolution of such systems on a three-dimensional Eulerian grid and provides an 
alternative to both N-body and phase space methods. 

Collisionless systems are ones in which the constituent particles move under the influence of 
the mean gravitational potential generated by all of the other particles. The state of a collisionless 
system is specified by a distribution function / = /(x, v,i) which gives the density of particles in 
phase space as a function of time. / is treated as a continuous fluid whose evolution is governed 
by the coupled Vlasov and Poisson equations (see e.g., Binney & Tremaine, 1986): 

W = y(Wdl_ Vi dl 
dt ~[ V dxi dvi 1 dxi 



E 

t=i 

V 2 K = AttG J d 3 vf . (2) 
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The Vlasov-Poisson pair must be solved numerically for most time-dependent systems of interest. 
One approach is to evolve / directly in phase space though efforts along these lines have had 
limited success primarily due to the large number of phase space dimensions typically involved. 
In contrast N-body or Particle methods (see, e.g., Hockney & Eastwood 1988) have been used 
successfully in virtually all astrophysical problems involving gravitational dynamics. The number 
of particles employed in an N-body simulation is typically many orders of magnitude less than 
the actual number of particles in the physical system one is modeling. This discrepancy can lead 
to unphysical effects due to "particle noise" such as two-body relaxation. Smoothing techniques 
alleviate this problem but at the cost of spatial resolution and care must be taken to choose 
a technique appropriate to the problem one is solving (Earn & Sellwood, 1995). On the other 
hand phase space methods, by explicitly constructing the smooth distribution function, avoid these 
difficulties. 

The SM attempts to circumvent the problems of both phase space and N-body methods by 
encoding phase space information in a continuous position space function which we call ijj. Given 
an initial distribution function at time U we can find the distribution function at a later time tf 
by the following general procedure: 

/(x, v, U) ip(x, U) Y>(x, t f ) /(x, v, t f ) . (3) 

Here DEM is the dynamical equation of motion for ip, M maps / to if), and M* maps ij) to /. 
Ideally M would be an invertable map with M* = M^ 1 in which case the procedure would yield 
an exact solution to the Vlasov-Poisson pair. 

With the SM, ip is formally identified as a complex Schrodinger field obeying the coupled 
Schrodinger and Poisson equations: 

ih^t = -—V 2 ^ + mV^ (4) 
ot 2m 



VV = 4vrGp . (5) 

M* is the coherent state or Husimi transform (Husimi 1940), essentially the square of a windowed 
Fourier transform, and M is a procedure to sample the initial distribution function to be discussed 
below. 

In N-body simulations fictitious "superparticles" are used to provide a statistical representation 
of the distribution function. Imagine dividing position space into cells smaller than the scales of 
interest but larger than the interparticle spacing. The number of particles in each cell gives the 
(course-grained) density field while their velocity distribution provides the remaining information 
contained in /. The resolution in position space relative to velocity space can be controlled by 
adjusting the cell size. However the size of a phase space resolution element, is set by the 
number of particles TV in the simulation: AQ ~ nCl/N. Here £1 is the volume in phase space filled 
by the system (6 dimensional in general; 3 dimensional for spherically symmetric systems with 
angular momentum (see below)) and n is the number of particles required to have a reasonably 
accurate estimate of / at a given phase space point. 

The SM shares some of these features. The de Broglie wavelength Ad e B for tp enters as a model 
parameter and is chosen to be as small as possible, in general a few grid spacings. Position space 
is divided into regions smaller than the scales of interest but larger than Ad c B so that the WKB 
approximation applies. We can then think of in this region as the superposition of plane waves 
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each of which represent different velocity streams in the distribution function. The correspondence 
principle guarantees that the "acceleration" of each stream is given by Newton's Law. Since ip is 
a continuous function there is no problem with particle noise. Moreover it is a function only of 
position rather than the phase space coordinates, and is therefore easier to work with. Formally, 
the resolution in phase space is set by the uncertainty principle. In practice it is set by the number 
of grid points N g : Af2 ~ n g £l/N g where n g is the number grid points across a typical de Broglie 
wavelength. We therefore expect the SM to be competitive with particle mesh simulations in terms 
of resolution in phase space as a function of the number of grid points. 

In a previous letter Widrow and Kaiser (1993) outlined the SM and presented simple illustrative 
simulations. In this paper we expand on that original work demonstrating, in a more quantitative 
way, the validity of the method. In particular we develop a new technique for handling "hot" 
systems (e.g. virialized systems) as well as systems with spherical symmetry. In Section 2 we 
review the SM and describe some of the new techniques implemented in this paper. Section 3 
presents results of test-bed simulations. In particular we simulate a stable equilibrium star cluster 
(Plummer sphere) and a nonequilibrium cluster. Both simulations assume spherical symmetry but 
include angular momentum. In Section 4 we make a few remarks concerning recent work on axion 
dynamics and in particular challenge claims made in the literature which suggest that the "field" 
nature of axions is important for understanding their behaviour in the Galaxy. Section 5 offers a 
summary and some concluding remarks. 

2 Schrodinger Method 

2.1 From ifj to JF 

We use the coherent-state (Husimi 1940) representation to construct a phase space distribution 
function .F(x, v) from t/j. Mathematically this is just the absolute square of a windowed Fourier 
transform. However, the physical content of this mapping is clearer when expressed in bra-ket 
notation: 

^(x,v) = |(r ? ( X ,v)|^)| 2 = |^(x,v)| 2 . (6) 
Here |r/(x, v)) is the wavepacket for a single 'particle' centered on the phase space point (x, v), 

/ 1 \ 3/2 / i \ 3/4 , „ 

<rKx,v)|x'> = (JLj (JLj e -(x-xo 2 /v-w.x 7 * (7) 

Equation (g) can then be seen as a prescription for decomposing ip into Gaussian wavepackets with 
velocity v and position x. One can show, directly from the Schrodinger equation, that T satisfies 

dT _A fSV dT dT\ 
dt ~{\ 9 x i dvi 1 dxi J 

(see, e.g., Skodje, Rohrs, & van Buskirk 1989) where L is the length scale over which the system 
varies and XdeB ~ IVV^^I — ft/ m v is the typically de Broglie wavelength of tp. This equation 
reduces approximately to the Vlasov equation provided AdcB *C rj <C L. 




+ o 



\ 2 

A deB 



~dt 



(8) 
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Moments of the distribution function are calculated in the usual way (e.g., (v 11 ) = J d 3 vv n F/ J cPvJ 7 ). 
In particular the density field p, which is the source term for the potential V, is given by 

p = jd\F= (^) 3/2 / d 3 x'e-( x " x ') 2 /" 2 |^(x')| 2 (9) 

Evidently p is found by applying a Gaussian filter to IV'I 2 ; hi effect, removing the high frequency 
modes. Of course what we are interested in is the potential V ~ V~ 2 /> and, since the inverse 
Laplacian also removes high frequency modes, we can safely replace p with in equation (||). 
In practice we calculate V from p using FFT. Particle mesh simulations also use FFT to calculate 
the force from p and so the CPU time for this part of the calculation will be the same for the two 
methods. 



2.2 From / to ip 

We begin by sampling the distribution function, compiling a list of phase space points much as we 
would in an N-body simulation. The initial ip is then taken to be the (incoherent) superposition of 
wavepackets centered on each of these points: 



1 N 



v El^ v i)> • ( 10 ) 

The \rj(xj,Vj)) are given by equation (^). The phase factor exp (iipj) is included to insure that the 
wavepackets will add incoherently. While the physical motivation for this prescription is clear (we 
are just summing up 'particles') a more mathematical motivation is the completeness of |ry(x, v)): 

J d 3 x J d 3 v|77(x, v))(r ? (x, v)| = / . (11) 

To see that this form provides a valid representation of / we compute ,F(x, v) = M* (Mf): 



-, N N 

•^( x > v ) = ^EE<«W (12) 

where 

Y. k = e -2imx-(vj-v fc )/R e -2imv-(xj-x fe )/a e 2i(m(x J -Vj-x fe -v fc )/?i+(^ J -v3 fe )/2) ^ 

Splitting this into a summation over equal and unequal indices we find, 

i N 

T = - e -( x -^) 2 /V e -2 m 2 ^ (v _ Vj) 2 M2 + ^ (i5) 

where T\ is a summation over the interference terms and will tend to be small. It is instructive to 
compare this result with the usual N-body distribution function: 

1 N 

/(x 1 v) = -^5(x-x J Mv-v J ). (16) 
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In the limit where the interference terms vanish equation (|T^) becomes equation (|jl|) with the delta 
functions smoothed into Gaussians. 

The discussion above reveals a key distinction between N-body methods and the SM: In an 
N-body simulation, smoothing is implemented only in the force calculation. The 'particles' in the 
SM are themselves treated as extended objects which evolve as would a distribution of particles 
(see below). Indeed, the particles/ wavepackets in the initial distribution mix with one another and 
do not keep their identity as do particles in an N-body simulation. We therefore have, in effect, a 
form of dynamical smoothing. 



2.3 Spherical Symmetry 

We now specialize to spherical symmetry where phase space has three dimensions; radial coordinate 
r, radial velocity v r , and angular momentum j = r^Jv 2 + v 2 ^. The Vlasov and Poisson equations 
take the form: 

df_ df_ , (f_<W\df_ _ () 
dt r dr I r 3 dr J dv r 

The derivation of equation ( |T7| ) assumes implicitly that j 7^ 0. (The Jacobian of the transformation 
from (x, v) to (r,v r ,j) is singular for j = 0). This in general does not present a problem so long 
as j is treated as a continuous variable. However for numerical work j is discretized and j = 
must be treated with care. Of course j can be viewed as a label since no differential operation with 
respect to j is ever performed. We can therefore write f(r, v r ,j) = fj(r, v r ) and treat the evolution 
of the different fj's separately. 



For radial orbits (j = 0) we write 



f (r,v r ) = ^llsMS^) . (19) 



where F satisfies the following equation: 



OF OF dV dF _ 

dt r dr dr dv r 



The Poisson equation is now 



d 2 dV . „ r°° 



dr dr 



' r ~~Q^ = ^Gj^r^F + n^fjj ■ (21) 



where the sum over j 2 is the discretized version of the integral over j 2 in equation (|i~g|). 

It is straightforward to extend the SM method to spherically symmetric systems with angular 
momentum. We construct a for each fj in the initial distribution function: For j = we have 



y (r,v r ) = C I dr'e-( r - r,) /2v e- imVrr ' /h Mr') (22) 
where tpo obeys the equation 



o 



di/j h 2 d 2 . <9Vo 



ih—— = ttV'o + rnVtpo 



dt 2m dr 2 dr 



= . (23) 

r=0 
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For j ^ 



*,-(r,v r ) = C r dr'e~ { - r - r 'y 2 l 2 ^e- imVrr 'l % (j) j {r') (24) 
o 



<9<f>- ?i 2 d 2 / 7 2 \ 

!r! ll = + m [ V + h) * ! *l- = ' (25) 

The distribution function is then 

r 2 \V \ 2 ~F ■ l^f-Zj (26) 
and the Poisson equation becomes 

-27T r2 7T V = i7T °P W 

where 

P = l^ol 2 + iEl^'l 2 - (28) 

P 

Equation ( p8| ) emphasizes an important point: the j = phase space map gives a density while 
the j ^ map produces dM(j 2 )/dr. 

2.4 Numerical Preliminaries 

It is convenient to write our equations in terms of the dimensionless quantities y = x/L, r = t/T, 
X = ^/\/7>0i an d U = j-zV where L and po are the characteristic size and density of the system of 
interest. We then have 

^ = -2^+-fX (29) 

vjc/ = |^^ 2 xx* (30) 

where a = hT/fiL 2 and 1 = 32GpoT 2 /3ir. We therefore have three dimensionless parameters, a, 
P, and 77/L, at our disposal. (3 determines the choice of timescale: T = j3Td where = y / 37r/32Gpo 
is the dynamical time. As we now show a determines the size of the phase space region accessible in 
the simulation while rj/L sets the relative resolution in position and velocity space. For simplicity 
consider one coordinate or two phase space dimensions. Suppose that our system comfortably fits 
inside a region of size C in units of L. We discretize position space by choosing yj = jC/N where 
N is the number of gridpoints. Further suppose that the system, when viewed in velocity space, 
comfortably fits inside a region of size V in units of L/T^. The lattice spacing in velocity space 
is connected to the lattice spacing in position space in the usual way (see, e.g., Press, Flannery, 
Teukolsky, & Vetterling 1986), i.e. the velocity is proportional to the wave number of the FFT. 
For discrete points in velocity space we have (in units of L/T) Uk = 2nka/£ which implies that we 
should set a = (3VC/2irN . Writing a = jC/N where 7 = O(l) we find (3 = 2-ir^/V. In practice we 
find that 7 ~ 5 works well. For the Plummer sphere simulations below we take V — 5 and j3 — 2n. 

The SM has the desirable feature that many quantities of interest can be calculated without 
having to actually construct the full phase space distribution function. We have already seen an 
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example of this when calculating p. In general, if the function does not depend on velocity, the 
velocity integral will collapse into a delta function. 

(Q) = |d 3 xd 3 vg(x) f(x,v) (31) 
= (2ir V 2 y 3 / 2 J d 3 x d 3 x' e -(x-x') 2 /^ 2 Q( x ) |^( x ')|2 (32) 
~ J d 3 x Q(x) |^(x)| 2 (33) 

In the case of velocity moments we can replace the v n term with a corresponding derivative, perform 
a change of coordinates and integrate by part to get, 



(34) 



s=0 



where is the unit vector in the j direction. For the important case of (v 2 ) this expression reduces 
to, 

{y2) ~ " (i)V d ' 3x {^ x ) V V(x)} , (35) 

where we have neglected a term proportional to ipip* , and have integrated by parts. This is just 
the usual quantum mechanical result (v 2 ) = (P 2 )/m 2 . 

Applying this to the case of spherical symmetry we then find, 

roo 

(Q) ~ ]T / dr Q(r) \^(r)\ 2 (36) 

One can, of course, resort to the 'brute force' calculation of more complicated quantities, 

(Q) =Y1 j dr J dVr Q( r i v r,3 2 ) Fj(r,v r )- (38) 
j 2 

The draw back of this approach is that it requires an 0(iV 4 ln(iV)) calculation for each j plane. 
These calculations can be very time consuming, and it is time well spent looking for a short cut of 
the kind just presented. 



3 Test-Bed Simulations 

3.1 Test Particles 

As a first concrete example of the SM we will consider the motion of 'test particles' in a fixed 
Plummer sphere potential, 

1/9 

V(r) = -GM (l 2 + r 2 ) . (39) 

The SM analogue of a test particle is a single coherent state wave packet centered at some chosen 
point in phase space. Here we will examine the motion of four particles, two in the j = plane 
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and two in the j = 0.5 plane. All four particles are taken to have zero initial radial velocity. 
The resulting initial wave functions are then just the sum of two Gaussians centered at the particle 
positions. The system is evolved by numerically solving equations (|j) and (p9|), using the algorithm 
presented in Visscher (1991). This is an explicit algorithm and is at least a factor of three faster 
than the usual implicit ones. Figures [I] and |2| show the four particles moving in the Plummer sphere 
potential. These plots help to show the physical connection between ip and J-, and the difference 
between the j ' = and j / phase space maps. The particle positions are given by the peaks 
of the wave function while the velocity information is encoded in the high frequency modes of ijj. 
From these plots one can see that the particles do follow the expected trajectories, justifying our 
choice of a spherical phase space map. 

These plots also allow an interesting comparison of the SM and N-Body methods at a funda- 
mental level. When one follows the evolution of an N-Body 'superparticle' it effectively carries with 
it a piece of phase space that has a fixed size and shape. When modeling systems in which mixing 
occurs these phase space blocks lead to unphysical coarseness of the distribution function. The 
evolution of the Schrodinger particles, on the other hand, carry with them a piece of phase space 
that has a variable shape and fixed volume. The motion of the test particles in figure |l] illustrates 
this point. The individual Schrodinger particles spread out along the phase space orbit in the same 
fashion as an equivalent distribution of point particles would. The conclusion is that fiducial SM 
particles can themselves take part in the phase mixing process. 

We can also examine the evolution of an initially 'cold' (i.e., single velocity stream) distribution 
of test particles. The corresponding wave function is given by square root of the desired density 
multiplied by an appropriate phase factor: 



V(x) = v / / 9(x)e-^( x ) (40) 
where V# = mv'. This leads to the phase space distribution 

■F(x,v) ~ p ^) e -v 2 (m V '-nveyn 2 _ ^ 

Compare this expression with the actual distribution function for a cold system: 

/(x,v)=p(x)5(v-v'(x)) : (42) 

The delta function has been replaced with a Gaussian of width ~ h/rj. 

Figure ^ shows the phase space evolution of initially cold (v r = 0) test-particles with angular 
momentum j = 0.5 (in units of L/Td) moving in the fixed potential of a Plummer sphere. The 
initial ip is given by 



dv r f(r,v r ,j 



2^ 



(43) 

3=0.5 



where / is the Plummer sphere distribution function (see below). With no velocity pressure to 
support itself the system quickly begins to phase-mix its way to equilibrium. 

3.2 Equilibrium Initial Conditions: The Plummer Sphere 

Our first full-scale test-bed calculation is of an equilibrium star cluster, specifically a Plummer 
sphere. The distribution function for this model is given by 

f{E) = { ^ if f <0 (44) 
Jy ' |0 otherwise v ' 
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where 

-1/2 

(45) 

is the energy per unit mass, A = ^^L 2 /G 5 M 5 , and M is the total mass of the system. 

Phase space is divided into Nj (r,v r ) planes of constant j with each plane having N 2 grid 
points. For the simulations that follow we set N = 1024 and Nj = 25, with no j = plane. 
(The distribution function is sharpley peaked in j 2 at the origin, and it is difficult to consistently 
integrate T with a small value for Nj . The lack of a j = plane in numerically consistent with a 
large value of Nj, where the j = plane would have a very small weight.) The planes are spaced 
equally in j although in general we can select planes at arbitrary values of j. We sample the set 
of Nj x N 2 , {r, v r }i phase space points using the usual N-Body technique (see e.g. Henon 1966). 
Figure |] illustrates this sampling method. Here the set of sampled phase space points are shown 
for the j = 0.5 and j = 1.0 planes. With these points in hand one can algorithmically compute 
the coherent state summation to produce the set of Nj initial wave functions. Figure ^ shows the 
results of this construction. Here a typical phase space plane is shown, along with the generating 
wave function and the mass and velocity distributions for that plane. One should note that with 
the above choice of 77, the r and v r resolutions are very similar. If r/ were increased, corresponding 
to a larger smoothing window, the mass distribution would become smoother while the velocity 
distribution would become more jagged. In terms of the Schrodinger particles of Figure |, the 
unit phase space particle would transform from a circle to a ellipse, becoming shrunken in the v r 
direction and elongated in the r direction. All of these effects would be transposed for r and v r if 
rj were, instead, to be decreased. 

The system is evolved by solving the coupled Schrodinger-Poisson system using the Visscher 
algorithm along with a potential solving routine based on Eastwood & Brownrigg (1979). Care must 
be taken to obtain the correct boundary condition at infinity when solving the Poisson equation 
on a finite grid. The routine that was implemented makes use of 'padding' to eliminate the effects 
of the grid boundary (at the expense of doubling the number of grid points). Figure |6| shows the 
mass and velocity distributions of the system in its initial state and after being evolved through 
10 dynamical times. Figure |7| shows the j plane of figure [| after 10 dynamical times. From these 
plots alone the system does appear to be stable. A quantitative measure of this is the virialization 
of the system: for a perfect Plummer sphere 2T + W = 0. Where T is the total kinetic energy: 

T=\{vl)+Y.f[(^)\n 2 dr, (46) 

and W is the potential energy of the system. The ratio 2T/\W\ was calculated during the evolution 
of the system over the course of 50 dynamical times. Figure I shows the evolution of 2T/\W\ 
for the system, with a similar curve obtained using a standard treecode (Barnes <fe Hut 1986) for 
reference. As can be seen from the plot, the SM Plummer sphere appears to be on par with the 
particle simulation. 



3.3 Non-equilibrium Initial Conditions 

As we saw in the last section, the Plummer sphere can be modeled quite well with only 25 j planes. 
But how well can this method model a perturbed system? One can imagine a situation where the 
numerical representation of a system is sufficient near equilibrium but is too coarse to accurately 
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model an evolving system. Therefore, in order to test the true dynamics of the spherical SM we 
performed a destabilized Plummer sphere simulation (Rasio, Shapiro, & Teukolsky (1989) use this 
system to test their phase space code). The system was destabilized by reducing the v r components 
of all the sampled Schrodinger particles by a factor of l/\/8, giving 2T/|W| = 17/24 ~ 0.708. This 
was also done for the BH simulation, again with iV = 1024 particles. Figure |9] shows the evolution 
of 2T/\W\ for both of these simulations. As can be seen, both undergo the expected bounce and 
then begin to virialize and settle into a new equilibrium. Figure 10 follows one of the j planes 
through some of this process to illustrate the phase mixing that occurs. 



4 Axion Dynamics and the Formation of Soliton Stars 

The SM is applicable to any collisionless system regardless of what form the constituent particles 
take. (The same holds for N-body and phase space methods.) This is particularly relevant for 
studies of galactic halos and large scale structure where the dominant component of the mass 
density, the so-called dark matter, is in some unknown form. Dark matter candidates range from 
astrophysical compact objects such as black holes and brown dwarfs to elementary particles such 
as neutralinos. Still the dynamics of the system is independent of nature's choice. 

There is at least one exception to the above argument. If dark matter is in the form of a 
very light scalar field, then "quantum-mechanical" effects will enter into the dynamical equations 
of motion. In this case, the Schrodinger-Poisson pair provides the exact equations of motion (for 
non-relativistic motion) for the system. The mass m is now a true physical parameter determined 
by particle physics theory. Of course so long as the Compton wavelength for the field is much less 
than the scales of interest in the problem, the field will behave like collisionless matter. 

Recently there has been considerable interest in soliton stars (for a review, see, e.g., Jetzer 
1992). These objects are essentially self-gravitating compact objects made up of bosonic fields. 
The interest in soliton stars has been generated largely by the conjecture that dark matter may be 
bosonic. Indeed one of the most popular dark matter candidates is the axion, a pseudo-Nambu- 
Goldstone boson that arises in the PQ (Peccei & Quinn 1977) solution to the strong CP problem 
(Weinberg 1978; Wilczek 1978). Axions that are viable (there are constraints on the axion mass 
from cosmology, stellar evolution, and supernova studies (see, e.g., Kolb & Turner 1990)) have a 
mass ~ 10~ 5 eV and therefore a de Broglie wavelength in a galactic halo of ~ 10 m. All this suggests 
that axions behave like any other dark matter candidate, i.e., like collisionless matter. However 
Seidel and Suen (1994) suggest that axions can form soliton stars. The question of whether this is 
indeed the case, or whether axions, like any other form of collisionless matter, form virialized clumps 
of matter, as in the axion miniclusters of Hogan and Rees (1988), is of fundamental importance in 
understanding whether or not axions are a viable dark matter candidate. Axions can annihilate 
into photons {AA — * 77) provided the density is high enough. In the axion miniclusters of Hogan 
and Rees (1988) the density is 10 g/cm and annihilation is unimportant. Alternatively the density 
in axion stars can be as high as 10 24 g/cm 3 . At these densities, the annihilation is very efficient 
making the configuration unstable. 

The question of course is whether axion stars ever form in the early Universe. A diffuse axion 
cloud must loose both angular momentum and energy if it is to form an axion star. Seidel and 
Suen (1994) find that a collapsing cloud can loose energy by radiating scalar material, a process 
they call gravitational cooling. However they only briefly mention the issue of angular momentum 
and carry out simulations of purely radial collapse. In these simulations the field collapses quickly 
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and an excited soliton star forms at the center and settles quickly by radiating scalar matter. 

In what follows we argue that any reasonable amount of angular momentum will prevent the 
object from reaching the stage where gravitational cooling can take place. Consider an initial cloud 
of mass M, radius Ri, binding energy E% ~ GM 2 /R, and angular momentum J. A useful quantity 
is the dimensionless parameter 

A - ^ (47) 

which is roughly the square root of the ratio between the rotational energy and the binding energy. 
As an object collapses, A oc R~ x l 2 provided there is no loss of angular momentum. We expect 
that this will indeed be the case at least in the initial collapse of the axion cloud. A/ < 1 in order 

1 /2 

for the final object to be gravitationally bound. We therefore require Xi(Ri/Rf) 1 < 1 where 
Rf — Ti 2 / n 2 GM is the radius of the final compact object. The bound is then 



where phaio = 0.008M Q /pc 3 is characteristic of the density of dark matter in the Galaxy. This 
is an incredibly small value for A. It is widely believed that galaxies acquire angular momentum 
during formation from tidal interactions with nearby mass distributions. Typical values from both 
analytic and numerical work suggest that A ~ 0.01 — 0.1 for galaxy-sized objects (see, e.g., Peebles 
1969; Efstathiou & Jones 1979). The conclusion is that axionic matter, if indeed it does make 
up the dark halo, will have too much angular momentum to ever reach the densities required for 
gravitational cooling and soliton star formation. 



5 Summary and Conclusion 

This paper describes test-bed simulations of collisionless, self-gravitating matter using the SM. 
Included are several improvements over the earlier work by Widrow and Kaiser (1993). In particular 
we describe how to set up "hot" (e.g., virialized) systems. We also develop the techniques required 
to handle spherically symmetric systems with angular momentum which may be useful for studying 
violent relaxation and equilibrium star clusters. Of course the method is not restricted to problems 
with high degrees of symmetry. Our focus on spherically symmetric systems is mainly for illustrative 
purposes and, in hindsight, may not have been the best choice because of the special features of 
spherical coordinates. Widrow and Kaiser (1993) ran 2D cosmological simulations with results that 
were in excellent agreement with those from N-body simulations. 

As discussed above, the SM is as efficient as a particle-mesh code in following the evolution of 
a system in phase space. One obvious improvement would be to use an adaptive grid and adaptive 
timestep to do the simulation. This could in principle make the method competitive with the 
widely used treecode of Barnes and Hut (1986). 

The SM can be adapted to a wide variety of problems that have, as their central equation, 
Vlasov. For example Widrow (1996) has shown that the SM (actually the Klein-Gordon equation) 
can be used to simulate collisionless systems in general relativity. The method can also be applied to 
problems in electrodynamics that require information about the phase space structure of particles 
thereby providing an alternative to particle- in-cell codes and Eulerian Vlasov codes (see, e.g., 
Manfredi et al. 1995). 
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N-body techniques have been used in virtually all areas in astrophysics that require numerical 
simulations of collisionless matter. However despite years of research, there are still questions that 
arise over the validity and applicability of these simulations. Most of the questions center around 
smoothing, the technique used to construct a smooth density field from a system of discreet particles 
(Earn & Sellwood 1995). The SM provides an alternative where the "super-particles" of the N-body 
experiment are replaced by a continuous field. While the field may look like the superposition of 
particle-like wavepackets (indeed, this is how we set up our initial conditions) these wavepackets 
spread out in phase space as the system evolves. We therefore have a type of dynamical smoothing. 
Whether this dynamical smoothing has any real advantages over N-body methods remains to be 
seen. 
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Figure Captions 



Fig. 1: Test Particles: Here the Schrodinger test particles are moving in a fixed Plummer 
sphere potential with j = 0. The solid cures are the expected phase space trajectories for the 
respective initial conditions. The top four plots are taken at T = and the bottom four at 
T = 10.2 Td- Here and below velocities are measured in units of L/Tj where Td = y / 37r/32G/>o (see 
text). 

Fig. 2: Test Particles: Here the Schrodinger test particles are moving in a fixed Plummer 
sphere potential with j = 0.5. The top four plots are taken at T = and the bottom four at 
T = 10.2 Td- The solid cures are the expected phase space trajectories for the respective initial 
conditions. 

Fig. 3: Cold Distribution: The initial phase space distribution corresponds to the j = 0.5 plane 
of a Plummer sphere with the velocity distribution suppressed. The system is out of equilibrium 
and rapidly collapses and begins to virialize. 

Fig. 4: Sampling the Plummer Sphere: This figure shows the j = 0.5 and j = 1.0 planes of the 
Plummer sphere , each with its 'N-Body' sampled particles. This sampling is the last step prior to 
the construction of the initial wave functions vpj. 

Fig. 5: The Equilibrium Plummer Sphere: This figure shows the initial wave function of a Plum- 
mer sphere j plane, along with the associated phase space distribution , Tj{r,v r ). The mass and 
velocity distributions were calculated directly from J-j. The solid curves are the model distributions. 

Fig. 6: The Equilibrium Plummer Sphere: This figure shows the initial mass and velocity 
distributions of the constructed Plummer sphere , along with the same distributions after the 
system was evolved through ten dynamical times. The dashed curves are the model distributions. 

Fig. 7: The Equilibrium Plummer Sphere: This is the j plane of figure || after the system was 
evolved through ten dynamical times. 

Fig. 8: The equilibrium Plummer Sphere: This figure shows the evolution of 2T/\W\ for the 
constructed Plummer sphere . The two curves are the results of the SM (line with points) and the 
BH tree code (with 1024 particles). 

Fig. 9: The Destabilized Plummer Sphere: This figure shows the evolution of 2T/\W\ for the 
destabilized Plummer sphere . The two curves are the results of the SM (line with points) and the 
BH tree code (with 1024 particles). 

Fig. 10: The Destabilized Plummer Sphere: This figure shows the evolution of j plane of the 
destabilized Plummer sphere . One can see that the system is phase mixing in a similar fashion to 
the cold system of figure ||. 
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